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Abstract 

We propose diffusion-like equations with time and space fractional derivatives 
of the distributed order for the kinetic description of anomalous diffusion and re- 
laxation phenomena, whose diffusion exponent varies with time and which, corre- 
spondingly, can not be viewed as self-affine random processes possessing a unique 
Hurst exponent. We prove the positivity of the solutions of the proposed equa- 
tions and establish the relation to the Continuous Time Random Walk theory. 
We show that the distributed order time fractional diffusion equation describes 
the sub-diffusion random process which is subordinated to the Wiener process 
and whose diffusion exponent diminishes in time {retarding sub- diffusion) leading 
to superslow diffusion, for which the square displacement grows logarithmically 
in time. We also demonstrate that the distributed order space fractional diffusion 
equation describes super-diffusion phenomena when the diffusion exponent grows 
in time {accelerating super- diffusion). 

PACS numbers: 

02.50.Ey Stochastic processes; 

05. 40. -a Fluctuation phenomena, random processes, noise and Brownian mo- 
tion. 



Recently, kinetic equations with fractional space and time derivatives have attracted 
attention as a possible tool for the description of anomalous diffusion and relaxation 
phenomena, see, e.g., [MKOO], [MLP01], [SZ97], [Mai96] and references on earlier stud- 
ies therein. It was also recognized [HA95], [Com96], [MRGSOO], [BMKOO] that the 
fractional kinetic equations may be viewed as "hydro dynamic" (that is, long-time and 
long-space) limits of the CTRW (Continuous Time Random Walk) theory [MW65] 
which was succesfully applied to the description of anomalous diffusion phenomena in 
many areas, e.g., turbulence [KBS87], disordered medium [BG90], intermittent chaotic 
systems [ZK93], etc. However, the kinetic equations have two advantages over a ran- 
dom walk approach: firstly, they allow one to explore various boundary conditions 
(e.g., reflecting and/or absorbing) and, secondly, to study diffusion and/or relaxation 
phenomena in external fields. Both possibilities are difficult to realize in the framework 
of CTRW. 

There are three types of fractional kinetic equations: the first one, describing Marko- 
vian processes, contains equations with fractional space or velocity derivative, the sec- 
ond one, describing non-Markovian processes, contains equations with fractional time 
derivative, and the third class, naturally, contains both fractional space and time deriva- 
tives, as well. However, all three types are suitable to describe time evolution of the 
probability density function (PDF) of a very narrow class of diffusion processes, which 
are characterized by a unique diffusion exponent showing time- dependence of the char- 
acteristic displacement (e.g., of the root mean square) [MKOO]. These processes are 
also called fractal, or self-affine processes, and they are characterized by the exponent 
H, called the Hurst exponent, which depends on the order of fractional derivative in 
the kinetic equation. We recall that the stochastic process x(t) is self-affine, or fractal, 
if its stationary increments possess the following property [ST94]: 

X(t + KT) - X(t) = K H (x(t + T) -X(t)) (1) 

where k and H are positive constants. The sign =implies, that the left and the right 
hand sides of Eq.(l) have the same PDFs. 

As a possible generalization of fractional kinetic equations, we propose fractional 
diffusion equations in which the fractional order derivatives are integrated with re- 
spect to the order of differentiation {distributed order fractional diffusion equations). 
They can serve as a paradigm for the kinetic description of the random processes pos- 
sessing non-unique diffusion exponent and hence, non-unique Hurst exponent. The 
processes with time-dependent Hurst exponent are believed to provide useful models 
for a host of continuous and non-stationary natural signals; they are also constructed 
explicitly [PV95], [AV99], [AVOO]. Ordinary differential equations with distributed or- 
der derivatives were proposed in the works by Caputo [Cap69], [Cap95] for generalizing 
stress-strain relation of unelastic media. In Refs. [BTOO], [BTOOa] the method of the 
solution was proposed which is based on generalized Taylor series representation. A 
basic framework for the numerical solution of distributed order differential equations 
was introduced in [DF01]. Very recently, Caputo [CapOl] proposed the generalization 
of the Fick's law using distributed order time derivative. 

We write the distributed order time fractional diffusion equation for the PDF f(x, t) 

as 
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/ = D % > /M) = *(*), (2) 





where r and D are positive constants, [r] = sec, [D] = cm2/sec, p(/3) is a dimensionless 
non- negative function of the order of the derivative, and the time fractional derivative 
of order f3 is understood in the Caputo sense [GM97]: 



df>f 



If we set p(/3) = 5(f3 — /3 ), < fi < 1, then we arrive at time fractional diffusion 
equation, whose solution is the PDF of the self-affine random process with the Hurst 
exponent equal to P /2. The PDF is expressed through Wright function [Mai97]. The 
diffusion process is then characterized by the mean square displacement 

oo 

(X 2 ) (t) EE J dxX 2 f( X ,t) = T{ p 2 o + 1) Drl ^ i,t(S °- ( 4 ) 
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This formula provides the generalization of the corresponding formula for classical 
diffusion valid at /3 = 1. Since (5 Q can be less than 1, Eq.(4) describes the process of 
slow diffusion, or sub- diffusion. 

Let us now prove that the solution of Eq.(2) is a PDF. The derivation here parallels 
to the method used in [SokOl]. Its aim is to show that the random process whose PDF 
obeys Eq.(2) is subordinated to the Wiener process. Returning to Eq.(2) and applying 
the transformations of Laplace and Fourier in succession, 



oo t 



f (k,s) = j dxe ikx J dte- st f(x,t) , (5) 



o 



we get from Eq.(2), 



where 
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I(sr) = J d(3{srfp{(3) (7) 
o 

We note that under the conditions described above the function I(st) is completely 
monotone on the positive real axis, i.e., it is positive and the signs of its derivatives 
alternate. We rewrite Eq.(6) as follows: 

oo oo 

/ (k, s ) = ~ J due- u ^ I+k2DT ] = J due' uk2Dr G(u, s) , (8) 



where 

G( u , s ) = l^p-e-" I W (9) 

is the Laplace transform of a function G(u, t) whose properties will be specified below. 
Now, f(x,t) can be written as 



l_XJ KJ\J 
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The function G(u, t) is the PDF providing the subordination transformation, from time 
scale t to time scale u. Indeed, at first we note that G(u,t) is normalized with respect 
to u for any t. With the use of Eq. (9) we get 
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(11) 



where L" 1 is an inverse Laplace transformation. Now, to prove the positivity of G(u, t), 
it is sufficient to show that its Laplace transform G(u, s) is completely monotone on 
the positive real axis [Fel71]. The last statement arises from the observation that 
G(u,s) is a product of two completely monotone functions, I/s and exp(-ul). The 
monotonicity of the former is obvious, whereas the monotonicity of the latter is an 
elementary consequence of the Criterion 2 in [Fel71], Chapter XIII, Section 4. Thus, 
we may conclude that the solution of Eq.(2) is a PDF, and that the random process, 
whose PDF obeys a distributed order time fractional diffusion equation, is subordinated 
to the Gaussian process using operational time. 

Since we are interested in diffusion problem, there is no need to obtain PDF, but 
its second moment only. Using Eq.(6) we get 



<- 2 > (t) 



d 2 f(k,t) 
dk 2 



k=0 



sI(st) 



Consider two fractional exponents in Eq.(2), namely, let 

p(P) = Brftf - + B 2 6(P - 2 ) , 



(12) 



(13) 



where < P 1 < P 2 < 1, Bl > 0, B2 > 0. Inserting Eq.(13) into Eq.(12) we get, 
denoting bi = B 1 T /3l ,b 2 = B 2 t^ : 



<x 2 > {t) = 2DTL; 1 
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Recalling the Laplace transform of the generalized Mittag-Lefffer function E^ u (z), fi > 
0, v > 0, which can be conveniently written as [GM97] 



L t {t»E^(-\t»} = 
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we get from Eq.(14), 



(15) 



(16) 



To get asymptotics at small t, we use an expansion, which is, in fact the definition of 
E^ u (z), see [Erd55], Ch.XVIII, Eq.(19): 



n=0 



r(fxn + v) ' 
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which yields in the main order for the square displacement , 
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For large t we use the following expansion valid on the real negative axis, see [Erd55], 
Ch.XVIII, Eq.(21), 

N -n 

W^-Ef^ + ^N" 1 -"). 1*1 (19) 

n=l x ' 

which yields 

^ K WftTT)(;J (20) 

Since /3 X < /3 2 , we have the effect of diffusion with retardation. We also note that 
the kinetic equation with two fractional derivatives of different orders appears quite 
naturally when describing subdiffusive motion in velocity fields [MKS98]. In this case 
the order of derivatives are (5 and (3 — 1, so that the situation differs from one discussed 
above. 

Now we consider a simple particular case which, in some sense, is opposite to the 
cases considered above. Namely, we put 

p(/3) = l,0</3<l. (21) 

Inserting Eq.(21) into Eq.(7) we get 

I(sr) = r — -, (22) 
log(sr) 

and, using then Eq.(12), 

(x 2 ) = 2Dr {log l - + 7 + e^ El (^j } , 
where 7 = 0.5772: is the Eiler constant, and 



(23) 



00 

y 



E 1 (z) = / dy— (24) 

j y 

z 

is the exponential integral. Using now the expansions valid on the positive real axis, 
see [AS65], Eqs.(5.1.11) and (5.1.51), respectively, 



El (z) = -7-logz-J2 { -^-, z^O, (25) 

n=l 

and 

-z N I 

V E^ 1 )"^' ( 26 ) 
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we get, retaining the main terms of the asymptotics at small and large times, respec- 
tively: 

^>~\ 2Drlog(±),f-oo (27) 

Thus, at small times we have slightly anomalous super-diffusion, whereas at large times 
we have superslow diffusion. 



The formula (27) can be generalized to the case 
P(0) = 



^- 1 ,o<p 1 <p<p 2 <i 

0, otherwise 



(28) 



and J dj3p(j3) = 1. Inserting Eq.(28) into Eq.(7) and, then into Eq.(12), we get 
o 
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Recalling Eq.(15) mean square displacement can be written as 
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Using Eq.(17), we get an expansion for (x 2 ) at small t: 
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where 



+V> (1 + & + n (/3 2 - &))} T- 1 (1 + /3 2 + n (/3 2 - &)) 



^(1/) = d (log !»)/*/ 



(30) 



(31) 



is the tp - function. At large t we explore the asymptotics valid on the real positive 
axis see [Erd55], Ch.XVIII, Eq.(22), 



E rAz) = V-V,^ (,./,) _ E ^— + O (|,|—) 

n=l v 

With using Eq.(32), Eq.(30) takes on the form 



2Dt ftV^fO~ n ^~ M 



{log - 



- V (1 + P 2 - n {(3 2 - P.mT- 1 (l + P 2 -n {(3 2 - &)) 



(32) 



(33) 



If we set /? x = 0, P 2 — 1 in Eqs.(31) and (33), then we arrive at the same expansions 
which are obtained by inserting Eqs.(25) and (26) into Eq.(23), respectively In par- 
ticular, the main term of the series (31) and (33) at (5 X — 0, (5 2 — 1 coincide with 
Eq.(27). 

The fractional diffusion equations with a given order of fractional time derivative 
are closely connected to the continuous-time random walk processes (CTRW) with 



the power-law distribution of waiting times between the subsequent steps [MKOO], 
[BMKOO]. Now we establish the connection between the distributed order time frac- 
tional diffusion equations and more general CTRW situations. Recall the basic formula 
of the CTRW in the Fourier-Laplace space [KBS87]: 

^f'* ' , (34) 
l-'p(k.s) 

A 

where w(s) is the Laplace transform of the waiting-time PDF w(t), and ip (k, s) is the 
Fourier-Laplace transform of the joint PDF of jumps and waiting times t). Assume 
the decoupled joint PDF, i/)(£,t) = \(£)w(t), and that the jump length variance is 
finite, that is, the Fourier transform of A (£) is 

X(k) « 1 - Drk 2 (35) 

to the lowest orders in k. Then, we consider the situations, in which mean waiting 
time diverges, that is, at large t the waiting time PDF behaves as 

w(t) « T p /t 1+ ^ < < 1, (36) 

and, consequently, 

w(s) « 1 - (srf (37) 

at small s. If f3 is constant, then inserting Eqs.(37) and (35) into Eq.(34) and making 
an inverse Fourier-Laplace transform, we arrive at time fractional diffusion equation. 
Now let us consider the case when (3 fluctuates. Indeed, for example, in the model called 
the Arrhenius cascade, which is inspired from studies of disordered systems, the unique 
P appears only under the assumption that the random trapping time is related to the 
random height of the well by the Arrhenius law [Bar99]. In a more realistic model, this 
law gives only the average value of the trapping time. Thus, we may speculate that in 
order to take into account the fluctuations of the trapping time, we can introduce the 
conditional PDF 

w(t\P) « T p /t 1+p , (38) 
and the PDF p{P), as well. Now, we have the relation 

i 

w (t) = [ dPp(P)w(t;P), (39) 



o 



where [0;1] is the whole interval for variations of p. We note that all waiting-time 
distributions with p > 1 correspond to similar behavior described by the first order 
derivative. Then, for the w(s) we have, instead of Eq.(37), 

i i 
w(s) w 1 - J dp(sTfp(P), p{p) > 0, J dpp(p) = 1. (40) 



Inserting Eqs.(40) and (35) into Eq.(34) we arrive at Eqs.(6) and (7). Thus, we see 
that the weight function p(P) has the meaning of the PDF. 



The model with fluctuating (3 is, of course, only one of the possible interpretations 
of Eq.(39): the non-exact power-law behavior of the waiting-time PDF can physically 
have very different reasons. In particular, the representation (39) allows us to consider 
regularly varying waiting- time PDFs, i.e., those which behave at t — > oo as w(t) oc 
t l f3 g(t) > where g(t) is a slowly varying function, e.g., any power of log t [Fel71]. 
We are also able to consider a waiting-time PDFs w(t) which show an approximately 
scaling behavior with the exponents changing with time. For such distributions the 
effective PDFs p((3) can be determined, and thus such non-perfectly scaling CTRWs 
can be described through distributed-order diffusion equations. The formal inversion 
of Eq.(39) can follow by noting that tw(t) taken as a function of log t is the Laplace- 
transform of the function 

JrmO<^<l 
nP) \ 0,K/?<oo ' 

Indeed, 

oo oo 

tw(t) = J dPcj>(P)t- p = J #W)exp( —P log t) — Lp {003) } (log t). 



The function 4>(f3) is thus given by 

m = K 1 {e u w(e u )}- 
The value of r can then be found through the normalization condition 



/ 



-p 



which defines then the function p((3). The description of the process through distributed- 
order diffusion equation is possible whenever this function is non-negative and concen- 
trated on < < 1. 

Now we turn to another type of fractional equation, namely, distributed order space 
fractional diffusion equation which, in dimensional variables, takes on the form 

2 

% = f daD{a)^L, f{x, 0) = 5(x), (41) 
at J d \x\ 

o+ 

where D is a (dimensional) function of the order of the derivative a, and the Riesz 
space fractional derivative d a /d \x\ a is understood through its Fourier transform $ as 

*G30 + -'* r/ - (42) 

If we set D(a) = K ao S(a — ao), then we arrive at the space fractional diffusion equation, 
whose solution is a Levy stable PDF of the self-affing stable process whose Hurst 
exponent is equal l/«o- The PDF is expressed in terms of the Fox's H-function [Fox61], 
[Sch86] . In the general case D (a) can be represented as 



D(a) = l a ~ 2 DA{a), 



(43) 



where / and D are dimensional positive constants, [/] = cm, [D] = cm 2 /sec, A is a 
dimensionless no n- negative function of a. The equation which follows for the charac- 
teristic function from Eq.(41) has the solution 



2 



/(M)=exp| ~ J daA{a){\k\l) a y (44) 
Note that the normalization condition, 

dxf(x,t) = f(k = 0,t) = l, (45) 



/ 



is fulfilled. 

Consider the simple particular case, 

A(a) = A^a - a ± ) + A 2 S(a - a 2 ), (46) 

where < a\ < a 2 , Al > 0, A2 > 0. Inserting Eq.(46) into Eq.(44) we have 

f(k,t) = exp{ -ai \k\ ai t-a 2 \k\ a2 t}, (47) 

where a x = A\D / l 2 ~ ai , a 2 = A 2 D / l 2 ~ a2 . The characteristic function (47) is the product 
of two characteristic functions of the Levy stable PDFs with the Levy indexes ai, a 2 , 
and the scale parameters a 1 /" 1 and a!/ a2 , respectively. Therefore, the inverse Fourier 
transformation of Eq.(58) gives the PDF which is the convolution of the two stable 
PDFs, 



oo 

f M = a -^W'-r^ J ^ (^) W 



(48) 



where L Q (x) is the PDF of the symmetric Levy stable law possessing the characteristic 
function 

4,o(A;) = exp(-|A;r). (49) 

The PDF given by Eq.(48) is, obviously, positive, as the convolution of two positive 
PDFs. The PDF will be also positive, if the function A(a) is represented as a sum of 
iV delta-functions multipolied by positive constants, iV is a positive integer. Moreover, 
if A(a) is a continuos positive function, then discretizing the integral in Eq.(41) by a 
Riemann sum and passing to the limit we can also conclude on the positivity of the 
PDF. 

Since the mean square displacement diverges for the Levy stable process, the anoma- 
lous super diffusion can be characterized by the typical displacement 5x of the diffusing 
particle [WS82], 

5xoc (\x\ q ) 1/q , (50) 

where (\x\ q ) is the q-th. absolute moment of the PDF obeying Eq.(41). For the stable 
process with the Levy index a 



x 



J C(q; a)t q / a , < q < a < 2 , s 

] oo, q > a 



where the coefficient 



C(q; a) = — (K a tf /a sin (^) T(l + q)T (l - i) (52) 
irq \ 2 / \ a J 

was obtained in [WS82]. To evaluate the g-th moment for the case given by Eq.(46), 
q < cti, we use the following expression , see, e.g., [Z0I86] 

00 

(|sH = ^r(l + g)sin(y) J dk(l-Ref(k,t) > )k- q - 1 . (53) 



We insert Eq.(47) into Eq.(53) and expand in series either exp(— a\ \k\ ai t), or exp(— a 2 \k 
with subsequent integration over k. As the result, for the g-th moment we have expan- 
sions valid at q < a\ and for small and large times, respectively, 

<|af ) = — (a 2 t)^ sin (^) T(l + q)F ( 1 - ±) x 
irq \ I / y a 2 / 




One can see, that at small times the characteristic displacement grows as t 1 /" 2 , whereas 
at large times it grows as t l / ai . Thus, we have superdiffusion with acceleration. 

In summary, we believe that the distributed order fractional diffusion equations can 
serve as a useful tool for the description of complicated diffusion processes, for which 
the diffusion exponent can change in the course of time. Further investigations are 
needed to establish the connection between proposed kinetics and multifractality. On 
the other hand, the development of numerical schemes for solving distributed order 
kinetic equations and for modelling sample paths of the random processes governed by 
these equations is also of importance 
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